R Notebook

setwd("~/Projects/BayesianStatistics/Chapter 16")

#16.1.4 Open your statistical software (R, Python, Matlab, and so on) and load any packages necessary to use Stan. (Hint: in R this is done by using library(rstan); in Python this is done using import pystan.)

library(rstan)

options(mc.cores = parallel::detectCores())

#16.1.5 Load the data into your software and then put it into a structure that can be passed to Stan. (Hint: in R create a list of the data; in Python create a dictionary where the ‘key’ for each variable is the desired variable name.)

discoveries <- read.csv('evaluation_discoveries.csv')
plot(discoveries$time,discoveries$discoveries, type='l', xlab = 'Year', ylab = '#discoveries')

X <- discoveries$discoveries
N <- length(Y)
dataList <- list(N=N, X=X)

#Problem 16.1.6 Run your model using Stan, with four chains, each with a sample size of 1000, and a warm-up of 500 samples. Set seed=1 to allow for reproducibility of your results. Store your result in an object called fit.

fit <- stan(file='discoveries.stan', data=dataList, iter=1000, chain=4, seed=1)

starting worker pid=73770 on localhost:11356 at 23:01:55.506
starting worker pid=73784 on localhost:11356 at 23:01:55.976
starting worker pid=73798 on localhost:11356 at 23:01:56.357
starting worker pid=73812 on localhost:11356 at 23:01:56.739

#16.1.7 Diagnose whether your model has converged by printing fit

print(fit)

Inference for Stan model: discoveries.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
lambda  3.10    0.01 0.18  2.74  2.98  3.10  3.22  3.46   817    1
lp__   39.84    0.02 0.73 37.61 39.66 40.13 40.31 40.36  1055    1

Samples were drawn using NUTS(diag_e) at Sun Aug 30 23:02:21 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

#16.1.9. Find the central posterior 80% credible interval for λ.

print(fit, pars='lambda', probs = c(0.1, 0.9))

Inference for Stan model: discoveries.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

       mean se_mean   sd  10%  90% n_eff Rhat
lambda  3.1    0.01 0.18 2.87 3.34   817    1

Samples were drawn using NUTS(diag_e) at Sun Aug 30 23:02:21 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls

lLambda <- extract(fit, 'lambda')[[1]]

Warning message:
In readChar(file, size, TRUE) : truncating string with embedded nuls

c(quantile(lLambda, 0.1), quantile(lLambda, 0.9))

     10%      90% 
2.874975 3.337161 

#16.1.10 Draw a histogram of your posterior samples for λ

params <- extract(fit)

hist(params$lambda)

#16.1.11 Load the evaluation_discoveries.csv data and graph it. What does this suggest about our model’s assumptions?

plot(discoveries$time, discoveries$discoveries, type='l')

hist(discoveries$discoveries)

#16.1.12. Create a generated quantities block in your Stan file, and use it to sample from the posterior predictive distribution. Then carry out appropriate posterior predictive checks to evaluate your model. (Hint: use the function poisson_rng to generate independent samples from your lambda).

fit <- stan(file='discoveries.stan', data=dataList, iter=1000, chain=4, seed=1)

DIAGNOSTIC(S) FROM PARSER:
Info: assignment operator <- deprecated in the Stan language; use = instead.

When you compile models, you are also contributing to development of the NEXT
Stan compiler. In this version of rstan, we compile your model as usual, but
also test our new compiler on your syntactically correct model. In this case,
the new compiler did not work like we hoped. By filing an issue at
https://github.com/stan-dev/stanc3/issues with your model
or a minimal example that shows this warning you will be contributing
valuable information to Stan and ensuring your models continue working. Thank you!
This message can be avoided by wrapping your function call inside suppressMessages()
 or by first calling rstan_options(javascript = FALSE).
Error in context_eval(join(src), private$context, serialize) : 
  0,248,Failure,-3,

starting worker pid=77647 on localhost:11356 at 14:01:45.915
starting worker pid=77661 on localhost:11356 at 14:01:46.211
starting worker pid=77675 on localhost:11356 at 14:01:46.503
starting worker pid=77689 on localhost:11356 at 14:01:46.771

SAMPLING FOR MODEL 'discoveries' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 1.5e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.012707 seconds (Warm-up)
Chain 1:                0.011052 seconds (Sampling)
Chain 1:                0.023759 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'discoveries' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 1.5e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.012615 seconds (Warm-up)
Chain 2:                0.011328 seconds (Sampling)
Chain 2:                0.023943 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'discoveries' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 1.5e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.012011 seconds (Warm-up)
Chain 3:                0.012002 seconds (Sampling)
Chain 3:                0.024013 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'discoveries' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 1.6e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.015401 seconds (Warm-up)
Chain 4:                0.011272 seconds (Sampling)
Chain 4:                0.026673 seconds (Total)
Chain 4: 

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls

lXSim <- extract(fit, 'XSim')[[1]]
lMax <- apply(lXSim, 1, max)
library(ggplot2)
qplot(lMax)

sum(lMax >= 12) / length(lMax)

[1] 0.0145

#16.1.13 A more robust sampling distribution is a negative binomial model:

        Xi ∼ NB(μ, κ) (16.2)
        

where μ is the mean number of discoveries per year, and var(X) = μ + μ2 . Here κ measures the κ degree of over-dispersion of your model; specifically if κ ↑ then over-dispersion↓.

fit_negbin <- stan(file='discoveries_negbin.stan', data=dataList, iter=1000, chain=4, seed=1)

starting worker pid=79491 on localhost:11356 at 14:39:36.746
starting worker pid=79505 on localhost:11356 at 14:39:37.096
starting worker pid=79519 on localhost:11356 at 14:39:37.424
starting worker pid=79533 on localhost:11356 at 14:39:37.796

SAMPLING FOR MODEL 'discoveries_negbin' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 3.1e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.077917 seconds (Warm-up)
Chain 1:                0.072575 seconds (Sampling)
Chain 1:                0.150492 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'discoveries_negbin' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 3.1e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.090217 seconds (Warm-up)
Chain 2:                0.071285 seconds (Sampling)
Chain 2:                0.161502 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'discoveries_negbin' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 3.7e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.37 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)

SAMPLING FOR MODEL 'discoveries_negbin' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 3.8e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.079128 seconds (Warm-up)
Chain 3:                0.059651 seconds (Sampling)
Chain 3:                0.138779 seconds (Total)
Chain 3: 
Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.076482 seconds (Warm-up)
Chain 4:                0.070567 seconds (Sampling)
Chain 4:                0.147049 seconds (Total)
Chain 4: 

There were 20 warnings (use warnings() to see them)

print(fit_negbin)

Inference for Stan model: discoveries_negbin.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

       mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
mu     3.12    0.01 0.22  2.72  2.97  3.11  3.26  3.59  1435    1
kappa  7.01    0.14 3.67  3.04  4.64  5.99  8.25 16.81   644    1
lp__  45.32    0.04 1.00 42.65 44.88 45.60 46.06 46.35   802    1

Samples were drawn using NUTS(diag_e) at Mon Aug 31 14:39:56 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

#16.1.14. Carry out posterior predictive checks on the new model. What do you conclude about the use of a negative binomial here versus the simpler Poisson?

fit_negbin <- stan(file='discoveries_negbin.stan', data=dataList, iter=1000, chain=4, seed=1)

DIAGNOSTIC(S) FROM PARSER:
Info: assignment operator <- deprecated in the Stan language; use = instead.

When you compile models, you are also contributing to development of the NEXT
Stan compiler. In this version of rstan, we compile your model as usual, but
also test our new compiler on your syntactically correct model. In this case,
the new compiler did not work like we hoped. By filing an issue at
https://github.com/stan-dev/stanc3/issues with your model
or a minimal example that shows this warning you will be contributing
valuable information to Stan and ensuring your models continue working. Thank you!
This message can be avoided by wrapping your function call inside suppressMessages()
 or by first calling rstan_options(javascript = FALSE).
Error in context_eval(join(src), private$context, serialize) : 
  0,248,Failure,-3,

starting worker pid=80068 on localhost:11356 at 14:45:10.832
starting worker pid=80082 on localhost:11356 at 14:45:11.136
starting worker pid=80096 on localhost:11356 at 14:45:11.404
starting worker pid=80110 on localhost:11356 at 14:45:11.677

SAMPLING FOR MODEL 'discoveries_negbin' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 5.8e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.127824 seconds (Warm-up)
Chain 1:                0.092388 seconds (Sampling)
Chain 1:                0.220212 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'discoveries_negbin' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 3.4e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.34 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.131062 seconds (Warm-up)
Chain 2:                0.108357 seconds (Sampling)
Chain 2:                0.239419 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'discoveries_negbin' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 3.8e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.116583 seconds (Warm-up)
Chain 3:                0.08936 seconds (Sampling)
Chain 3:                0.205943 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'discoveries_negbin' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 3.4e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.34 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.100734 seconds (Warm-up)
Chain 4:                0.084441 seconds (Sampling)
Chain 4:                0.185175 seconds (Total)
Chain 4: 

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls
9: In readChar(file, size, TRUE) : truncating string with embedded nuls
10: In readChar(file, size, TRUE) : truncating string with embedded nuls

lXSim <- extract(fit_negbin, 'XSim')[[1]]
lMax <- apply(lXSim, 1, max)
library(ggplot2)
qplot(lMax)

sum(lMax >= 12) / length(lMax)

[1] 0.2515

16.1.15

Find the central posterior 80% credible interval for the mean rate of discoveries μ from the negative binomial model. How does it compare with your results from the Poisson model? Why is this the case?

print(fit_negbin, pars='mu', probs = c(0.1, 0.9))

Inference for Stan model: discoveries_negbin.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

   mean se_mean   sd  10%  90% n_eff Rhat
mu 3.13    0.01 0.22 2.84 3.42  1334    1

Samples were drawn using NUTS(diag_e) at Mon Aug 31 14:45:49 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

#16.1.16. Calculatetheautocorrelationintheresidualsbetweentheactualandsimulated data series. What do these suggest about our current model?

lXSim <- extract(fit_negbin, 'XSim')[[1]]
mResid <- sweep(lXSim, 2, X)
lCorr <- sapply(seq(1, 200, 1), function(i) acf(mResid[i, ], lag.max=1)$acf[[2]])

qplot(lCorr)